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Abstract 

We investigate the thermal expansion of crystalline SiC>2 in the j3- cristo- 
balite and the /3-quartz structure with path integral Monte Carlo (PIMC) 
techniques. This simulation method allows to treat low-temperature quan- 
tum effects properly. At temperatures below the Debye temperature, thermal 
properties obtained with PIMC agree better with experimental results than 
those obtained with classical Monte Carlo methods. 

I. INTRODUCTION 

The study and analysis of materials properties of crystalline silicates are important since 
these systems are used in many industrial processes and they occur also in many natural 
rocks. Many interesting effects have been found by experimental techniques at temperatures 
well below the Debye temperature. In this temperature range quantum effects like zero point 
motions and corresponding derealizations of atoms are important which have to be taken 
into account in serious theoretical studies. Usually this is done by lattice dynamics theories, 
in the framework of the harmonic or quasiharmonic approximation. However, both near 
second-order structural phase transitions and quite generally at higher temperatures the 
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accuracy of this approach is sometimes uncertainly, and methods that work at all conditions 
would be desirable. 

PIMC simulationsi enable us to analyse the crystal low temperature thermal properties. 
In principle, this method yields exact quantum-statistical averages (apart from statistical 
errors) and reduces to classical statistical averages at high temperatures. In general the 
agreement with experimental data is better compared to classical computations. For ft- 
cristobalite it turns out that even at temperatures as high as 600 K only with PIMC a 
good agreement with experimental findings is obtained. The negative thermal expansion in 
/5-quartz however may be explained by a classical modelling. 

The outline of this paper now is as follows: in Sec. II we define the models that are 
studied, in Sec. Ill we briefly review the PIMC technique, while Sec. IV describes our 
results in detail. We conclude in Sec. V with a brief summary. 



II. THE MODEL 

Starting from the diamond lattice of silicon solids and placing between each neighbour 
pair of Si atoms an O atom, we arrive at a silicate which exists in nature in cubic symmetry: 
ft- cristobalite. 

For the computation of Si02-structures a variety of potentials are studied in the litera- 
ture, for a review see reffi Concerning the computational time it is of great advantage that 
apparently two-body potentials like the TTAM-potential@ and the BKS-potentiali describe 
the system properties at least as good as three body potentials. Thus most of the investiga- 
tions presented here have been done with the TTAM- and the BKS- potentials. Compara- 
tive studies have been done as well with other potentials like the two-body contribution of 
a complicated potential obtained via ab-initio computations^ ("core-shell" -potential). 

In case of the TTAM and the BKS potentials the interactions between the Si and O- 
ions are modelled by pair potentials of the form: 
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$( r . .) = + A .. eX p(-^r lJ ) - ^ (1) 

r »j r ij 

Here e 2 = 4 1 8 6 g 2 4 1 2 9 7r e\^A. In Eq.(|I|) each ion carries a partial charge: A Si- ion has the charge 
Q = 2.4 and an O- ion has the charge Q = —1.2. 

The TTAM potential has been developed by Tsuneyuki et.alfi The parameters in Eq.(|l]) 
are listed in table I. The parameters in Eq.(|l|) for the BKS- potential! are tabulated in 
table II. 

In general complex materials with several components interacting with long ranged 
Coulomb interactions have to be handled by appropriate numerical procedures such as the 
Ewald summation technique. In the numerical simulation the computation of the Coulomb 
contribution to the energy requires a summation not only over the particle indices but also 
over all periodic boxes: 

(iy^j for n— 0) 

The sum in Eq.(|D is over all periodic images shifted by n and all particle indices except 
i — j for n = 0. For numerical details of the computation of the Coulomb energy and how 
to use a relatively simple and efficient algorithm we refer to Ref.@ 

Both potentials have in common that at small interparticle distances the atomic repulsion 
is neglected and rather a divergency in the potentials at small distances appears. In order 
to avoid a "Coulomb-collapse" of the system, at small distances we have to include an 
additional repulsive term. This term can be modelled as a harmonic contribution! suitable 
for molecular dynamics simulations or as a purely hard sphere repulsion. In this study for 
the sake of simplicity we follow the latter prescription, the hard sphere distance is chosen as 
the distance of the first potential maximum in Eq. (0). 

III. PATH INTEGRAL MONTE CARLO SIMULATIONS 



Canonical averages < A > of an observable A in a system defined by the Hamiltonian 
Ti = Ekin + Vp t of iV particles in a volume V are given by: 
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<A>=Z~ X Sp [Aexp(-pH)} . (3) 

Here Z = Sp [exp(— /3TC)] is the partition function and (3 = l/k B T is the inverse tempera- 
ture. Utilising the Trotter-product formula, 

exp((3H) = lim (exp(-/3E kin /P)exp(-/3V pot /P)) p , (4) 

P^oo 

we obtain the path integral expression! for the partition function: 

Z(N, V,T) = ^j£^) 3NP/2 / 4rW}exp[-^(f:^(r? - r^) 2 + ^({rW}))] 

(5) 

(s) 

Here, m is the particle mass, integer P is the Trotter number and r k denotes the coordinate 
of particle k at Trotter-index s, and periodic boundary conditions apply, the particle with 
Trotter-index P + 1 is the same as the particle with Trotter-index 1. This formulation of 
the partition function allows us to perform Monte Carlo simulations for increasing values of 
P approaching the true quantum limit for P — > oo. Note that Eq. (Q) does not take into 
account any quantum- mechanical exchange between particles (for atoms as heavy as Si and 
O this is an excellent approximation, though it would not work for He crystals). 

Thermal averages in the ensemble with constant pressure p are given via the correspond- 
ing partition function 

A(N,p,T)= dVexp[-PpV]Z(N,V,T) (6) 

<j 

. In order to make our results comparable with those of experiments, our PIMC simulations 
have been performed in the constant pressure ensemble. Most of our simulations! have been 
done with N = 684 particles (/3-cristobalite) , iV = 576 (/?-quartz), p = 0, and P values 
up to 100. A typical PIMC data point in figure || required about 1200 CPU hours on a 
CRAY-T3E (single processor). 

IV. RESULTS 
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A. Finite size effects and pair potential sensitivity 



In order to analyse effects due to the finite system size simulations for the system in 
the (3- cristobalite structure have been done using the three potentials described above and 
N = 648 particles as well as with the BKS- potential and iV = 192 particles. In Figs. [1] and 
|2] the potential energy and the systems box volume is shown. 

The energy obtained by the core-shell potential (in the two-body reduced form) is signif- 
icantly smaller than the energy obtained by the other two potentials. The volume obtained 
by the core-shell potential is much larger than the experimental data0 and the thermal 
expansion coefficient is negative over the entire temperature range studied- in contrast to 
the experimental results for /3-cristobalite. 

For the BKS- potentials the successful description of Si0 2 -melt properties! using a cut- 
off distance of 5.5A for the short ranged interactions suggested the choice of this cutoff- 
distance in the present study as well. From figure p] however it is apparent that the volume 
decreases strongly with increasing temperature- in contrast to the experimental data. For 
the ordered crystal it turned out that a larger cutoff value (12 A) is more adequate for the 
description of the system properties. 

The energies obtained by the TTAM- potential and by the BKS- potential (for both sys- 
tem sizes) agree well. In case of the volume however the results for the BKS- potential agree 
for both system sizes within numerical scatter and show a qualitatively similar temperature 
dependence as the volume obtained by the TTAM- potential. The quantitative values of 
the volume for the TTAM- potential however are much larger than the experimental data, 
whereas the volume obtained by the BKS- potential agrees well with the experimental data 
for temperatures above 700 K. 

In summary the best agreement with the experimental data can be obtained with the 
BKS- potential - at least in a temperature range above 700 K. For smaller temperatures 
apparently none of the classical computations can successfully describe the experimental 
data. 
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B. Quantum simulations for (3- cristobalite 



1. Potential energy and volume 



A different picture emerges if results of quantum simulations are considered in addition 
(BKS- potential, N = 648), see figure [|. At temperatures below 700 K the quantum results 
for the potential energies deviate from the classical results to higher values, but for higher 
temperatures they agree with the classical values. The volume of the simulation box obtained 
by the PIMC simulations (Fig. 0) agree within numerical scatter with the experimental data 
for all temperatures studied- in contrast to the classical simulations. Despite the apparently 
high temperatures of 700 K quantum effects are very important for the thermal properties 
of /3-cristobalite and thus quantum effects should not be neglected in simulations of real 
materials. This supports the assumption that quantum effects may still play an important 
role above the Debye temperature which is at about 500 K for SiC>2. 



middle of a bond connecting two Si- atoms, the Si-O-Si bond angle is 180° and the Si-0 



which of these models describes the properties of (3- cristobalite best is not yet decided by 
experimental methods. 

In figure ^| we show the bond length distribution function for the different atom pairs for 
three temperatures using the BKS potential and in figure f| the average bond lengths as a 
function of temperature. Note that due to the choice of a finite size cubic simulation box 
with periodic boundary conditions, /3-cristobalite is stable (or metastable, respectively) in 
the simulation over the entire temperature range, unlike experiment. 



2. Structural details 



According to a model of Wyckoflffl for (3- cristobalite the O- atoms are located in the 
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In the PIMC simulations (with P=30) the bond lengths have larger values at low tem- 
peratures compared to the classical case (P=l) (Fig. f|). The quantum results for the Si-0 
bond length are in much closer agreement with the literature value of 1.61 A. The Si-0 as 
well as the 0-0 bond lengths scale approximately linearly with the temperature, the Si-Si 
bond length decrease is stronger in parallel to the temperature dependence of the volume. 
This may be due to the fact that the positions of the Si-atoms fix the lattice constant of (3- 
cristobalite and the O- atoms have a much larger mobility in the crystal lattice. In figure |] 
we find again the different behaviour of the Si-Si bonds compared to the other bonds. The 
temperature dependence of the Si-0 and the 0-0 bonds have the following features: 

• The width of the distributions increases with temperature, the height decreases. 

• The average values increase with the temperature. 

• At the same temperature the distributions obtained by PIMC (P=30) have larger 
average values than the classical distributions. 

• The difference in the distributions and the average values between the classical and 
the quantum results decrease with the temperature. 

All these effects are present in the Si-Si bonds as well, but there is a significant difference: 
In Figs. ^| and ^ we find that for the other two bond types the differences between the classical 
and the quantum values at the temperature T are always smaller than the difference between 
the classical values at T and the next higher temperature studied. For the Si-Si bond the 
quantum distributions at 400 K are closer to the distributions at 1600 K (which do not 
differ very much any more) than to the classical distributions at 400 K. This shows that in 
particular the Si-Si bond properties are significantly determined by quantum effects. 

In order to obtain more information on the long-range molecular structure, the radial 
distribution functions g(r) were computed up to a distance of 10 A for Si-Si, Si-O, and 
0-0 bonds. The knowledge of g(r) allows the calculation of the so-called total correlation 
function T(r), which is a superposition of the various ^(r)'s@ T(r) has been measured 
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recently by neutron diffract ion. Ila The measurements allowed to accurately estimate the real 
Si-0 bond lengths and it was concluded that the typical Si-0 bond length is considerably 
longer than the Si-0 bond length in the so-called ideal /3-cristobalite geometry. In Fig. (|5]a), a 
comparison between the experimentally observed and the calculated function T(r) is shown. 
The agreement between the two curves is at least semi quantitative and the conclusions 
drawn in Ref.0 could be confirmed with the exception of the width of the first Si-0 peak and 
the first 0-0 peak. These peaks are broadened in experiment with respect to simulations, 
which include quantum effects. The broadenings might be due to lattice defects or to internal 
surfaces in the experimental samples. Both kind of defects were absent in our simulations. 
Similar results on T(r) were reported previously by purely classical simulationsll! on the 
basis of the TTAM potential. 

Fig. (||b) allows to relate the features in T(r) to the various g(r)'s. It is interesting to 
note that the peak in T(r) at about 5 A is due to a simultaneous maximum in goo{ r ) and 
5 , Sio( r ) ; while the peak in T(r) at 6.25 A does not have a corresponding peak in any g(r). 
The large local maximum in T(r) at about 9 A is then located at a position where all g(r)'s 
have a local maximum as well. The shape of T(r) is particularly sensitive to details of the 
potentials at large distances, e.g., cutting off the short-range part of the potential at 4.5 A 
alters T(r) significantly for r > 8 A. 

C. Quantum simulations for /3-quartz 

The (3- quartz structure differs from the (3- cristobalite structure by the positions of the 
atoms in the unit cell. (3- cristobalite has a cubic structure, (3- quartz has a hexagonal 
structure and at a temperatures below 846 K a trigonal a- quartz is stablSHIHl. For 
structural details of the (3- quartz see appendix A. The piezoelectricity (and pyroelectricity) 
of /3-quartz is at least reduced for (3- cristobalite due to symmetry reasons. 

As for (3- cristobalite the simulations have been done for Trotter values of P=l and 
P=30 at temperatures between 600 K and 1600 K. Generalizing Eqs. (3), (5) and (6), in our 
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constant pressure simulations a hydrostatic pressure p has been applied and fluctuation trials 
of the box lengths in x-, y- and z- directions have been attempted independently keeping 
the angles between the box axes constant. Lattice constants (a and c) have been obtained 
from the resulting average box geometry. The experimentally observed transition at 846 K 
is suppressed in the simulations since we keep our box in rectangular shape. A trigonal 
a- quartz structure thus cannot be stabilised, rather a metastable, undercooled /3-quartz is 
present in the simulations. Above the phase transition temperature, the structure in the 
simulation agrees with the structure deduced from the experimental data. 

In figure |6| the potential energy and lattice constants of /3-quartz are compared to ex- 
perimental data§. The energy has a linear temperature dependence in the classical case 
and slightly higher values in the quantum simulations, the difference to the classical values 
increases the lower the temperature. Apparently in this temperature region quantum effects 
become important- as can be expected due to a Debye- temperature of 500 K for SiC>2. The 
temperature dependence of the volume only qualitatively agrees with the experimental data. 
In figure |6| we note the qualitatively different behaviour between the experiments and the 
simulation data at temperatures close to the a ^ (3 transition at 846 K. As explained above 
this different behaviour is caused by the box shape constraint in the simulation. 

As in the experiment the lattice constant in a- direction is approximately temperature 
independent at high temperatures whereas in c- direction a negative thermal expansion coef- 
ficient is found. In the case of the lattice constant c, the PIMC results are in somewhat better 
agreement with the experimental data than the classical results. Since for /3-cristobalite it 
turned out that by using the literature value (5.5 A) for the cutoff the agreement with the 
experimental data is not as good as without the cutoff (see above), we chose the maximum 
cutoff value which is possible in our box geometry. The deviation of 1% (c-direction) and 
0.3% (a- direction) from the experimental values may be due to the choice of parameters 
in the interaction potential, in particular to the cutoff of the short ranged interactions in 
the BKS- potential. Probably a restriction to smaller cutoff values would result in a better 
agreement with experimental data. 
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A typical simulation data point in figure |] also required a computational effort of about 
1200 CPU hours on a CRAY-T3E computer (in units of single-processor time). 

Surprisingly the c- lattice constant obtained by a classical simulation decreases slightly 
with increasing temperature. This behaviour- which is present in the experimental data as 
well- may thus be explained by a purely classical argument. In the literature the 'rigid unit 
mo des&H are suggested as being responsible for the negative thermal expansion coefficient. 
For the sake of a more complete discussion of our results, we recall this argument here: 
The idea is that in a covalently bonded crystal structure one can identify units which are 
essentially rigid, because the length of the covalent bonds are essentially rigid, but that at 
certain atomic positions where the rigid units are linked together, they may rotate relative to 
each other. These modes may be visualised schematically as in figure [j] for two dimensional 
structures. The quadratic bases in this example may correspond to SiC>4 tetrahedron in 
(3- quartz. If the angle stabilising forces in such units are significantly stronger than the 
forces at the common O- atoms connecting such tetrahedron, thermal excitation leads to 
a rotation of the complete bases. In contrast to the positive thermal expansion coefficient 
caused by the anharmonicities of the interatomic interactions such collective modes may 
cause a negative thermal expansion coefficient. The possible mechanism is briefly sketched: 
By tilting of the unit cell by the angle the system volume A = A (see Fig. |7|) changes 
to A = t4 cos 2 (G). For small angles and harmonic rotational oscillations with frequency uj 
of the units (i.e. tetrahedra) with moment of inertia / the resulting thermal average of A 



The thermal average of the tilt angle < O 2 >t increases with temperature and thus leads 
to a decrease of system volume. At low temperatures this effect may more than compensate 
the 'normal' thermal expansion caused by the motion of single atoms such that the thermal 
expansion coefficient is negative. 
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1. Structural details (3 -quartz 

In order to get more information on the molecular structure, the radial distribution 
functions g(r) were computed for /3-quartz, again up to a distance of 10 A. The correlation 
function T(r), which follows from the various g(r)'s, is shown in Fig. (|8]a) together with 
9SiSh 9SiO, and in Fig. (|b). 

It is instructive to compare Fig. (H) for /3-quartz with the corresponding figure for (3- 
cristobalite, Fig. (^). The location of the peaks in g(r) are nearly identical for nearest 
neighbor Si-Si, Si-O, and 0-0 peaks, only the peaks are broadened in the case of /3-quartz, 
which can entirely be attributed to the larger temperature of the /3-quartz simulation. It 
can be concluded that the tetrahedra formed by Si04 units are basically identical for both 
lattice structures. While the radial distribution functions differ markedly between the two 
lattices for distances r > 3.5 A, there is a much less pronounced, yet noticeable, effect in 
T(r). E.g., T(r) shows a shoulder for /3-quartz at r ps 3.7 A which is absent in the simulation 
of /3-cristobalite. 

D. "Size" of the quantum particles 

By determining the average radius of the quantum-particles across the quantum chain 
their "size" can be mapped as a function of temperature. In figure [5] we compare the 
simulation results for these radii in the various lattice structures with the thermal de-Broglie 
wave lengths for free O- and Si- atoms. Apparently this quantity does not depend much 
on the particular lattice structure. This may be due to the fact that the quantum chain 
radius is determined by the close surrounding of the particles and the potentials. In both 
crystalline structures studied here the basic unit is a Si04 tetrahedron and the BKS potential 
is used. The de-Broglie wave lengths for the free particles are much larger than the radii of 
the particles of the crystals. This may be due to the stronger binding and localisation of 
the latter. 
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V. SUMMARY 



In this paper we present the results of our PIMC simulations for crystalline silicates. With 
PIMC simulations we are able to analyse the thermal properties of crystals at temperatures 
at which quantum effects have to be considered as well as anharmonic effects. In general 
the agreement with experimental data is better compared to classical computations. For 
/3-cristobalite it turns out that even at temperatures as high as 600 K only with PIMC a 
good agreement with experimental findings is obtained. The negative thermal expansion in 
/5-quartz however may be explained by classical effects. 

We conclude that for a realistic materials modelling the simulation by path integral 
Monte Carlo techniques is a useful method which should be used in future studies of ther- 
mal properties of other crystals as well. 
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VI. APPENDIX A 



In most simulations a rectangular simulation box is used. It is thus of advantage to map 
the hexagonal system onto a larger rectangular box such that the periodic order results in 
the hexagonal structure.! 



The lattice structure for (3- quartz is sketched in Fig. 10 and may be described as follows. 
In the space group P6222 the unit vectors are: 
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ei = (4.9977, 0, 0) 

e 2 = (-2.49885,4.328135,0) 

e 3 = (0,0,5.4601) 

(8) 

Thus the lattice constants are: 

a = 4.9977 A 
b = V3a 
c = 5.4601A 

(9) 

The coordinates of the particles in a unit cell are given in table III, the Cartesian coordinates 
in table IV. 

Due to the differences in box lengths one has to reduce the effect of different surface 
areas on the simulation results. This is done by matching sufficiently many unit cells such 
that the simulation box is approximately cubic. 
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FIGURES 
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FIG. 1. Classical (P=l) unit cell volume for /3-cristobalite. The simulation box is a cube 
containing 27 unit cells with N = 648 particles, except otherwise noted. Comparison with experi- 
mental data taken from Ref.0. 
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FIG. 2. Potential energy (a) and unit cell volume (b) for /3-cristobalite, results of path 
integral Monte Carlo simulations (with N = 648 particles and 27 unit cells unless otherwise noted). 
Comparison with experimental data taken from Ref.0. In addition simulations have been done 
with a "core-shell" potential! (see text) with an average potential energy of about -43.2 eV. 
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FIG. 3. Distributions of the bond lengths in (3- cristobalite at the temperatures T = 400 K 
and 1600 K. (a) Si-Si bonds, (b) Si-0 bonds, (c) 0-0 bonds. Comparison of classical (P=l) and 
quantum simulations (P=30). 
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FIG. 4. Bond lengths between nearest neighbours in (3- cristobalite. (a) Si-Si bonds, (b) Si-0 
bonds, (c) 0-0 bonds. Comparison of classical (P=l) and quantum simulations (P=30). 



19 




FIG. 5. a) T(r) for /3-cristobalite as a function of distance r at 573 K. Experiment is represented 
by a solid line and simulation data is represented by a dashed line, b) Corresponding radial 
distribution function g(r) for Si-Si, Si-O, and 0-0 bonds. The curves are normalized such that 
g{r) — ► 1 for r — ► oo. 
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FIG. 7. Schematic picture of rigid unit modes in two spatial dimensions. By rotation of the 
unit cells (squares) by the angle 6 the volume of the system can be reduced. Comparison of the 
system volume Aq with = (left side, dotted line) with the system volume A with 6/0 (right 
side, dashed line). 
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r[A] 



FIG. 8. a) T(r) for /3-quartz as a function of distance r at 900 K. b) Corresponding radial 
distribution function g(r) for Si-Si, Si-O, and 0-0 bonds. The curves are normalized such that 
g(r) — ► 1 for r — > oo. 
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FIG. 9. Quantum chain radii in (5- cristobalite and j3- quartz. Comparison with the de- Broglie 
wave lengths of free particles 
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TABLES 



ion pair 


A/eV 


B/A' 1 


C/eVl 6 


Si-Si 


872.4 


15.22 


23.265 


0-0 


1756.98 


2.846 


214.75 


Si-0 


10722.23 


4.796 


70.739 



TABLE I. Parameters in Eq.(|l|) for the TTAM-potential. 



ion pair 


A/eV 


B/A- 1 


C/eVl 6 


Si-Si 










0-0 


1388.7730 


2.76000 


175.0000 


Si-0 


18003.7572 


4.87318 


133.5381 



TABLE II. Parameters in Eq.(jT]) for the BKS-potential. 
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atom ' 


X 


y 


z 


Si 1 


1/2 





o 


Si 2 





1/2 


2/3 


Si 3 


1/2 


1/2 


1/3 


1 


X 


2x 


1/2 


2 


X 


2x 


1/2 


3 


2x 


X 


1/6 


4 


X 


X 


5/6 


5 


X 


X 


5/6 


6 


2x 


X 


1/6 



TABLE III. Coordinates of the particles in the unit cell of /3-quartz (see Fig. 10). 



atom 


x/A 


y/A 


z/A 


Si 1 


2.498850 








Si 2 


-1.249425 


2.164068 


3.640067 


Si 3 


1.249425 


2.164068 


1.820033 


1 





1.793579 


2.730050 


2 


2.498850 


2.534556 


2.730050 


3 


0.945565 


3.431346 


0.910017 


4 


-.945565 


3.431346 


4.550083 


5 


3.444415 


0.896790 


4.550083 


6 


1.553285 


0.896790 


0.910017 



TABLE IV. Cartesian coordinates of the particles in the unit cell of /3-quartz. 
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